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Abstract 

^ha  purpose  of  this  report  is  to  introduce  an  adaptive 
estimation  and  parameter  identification  scheme  »rhich  we  eheli  call 
Multiple  ^del  Estimation  Algorithm  (MMBA) . The  MMEA  consists 
of  a bank  of  Kalman  filters  with  each  matched  to  a possible 
parameter  vector.  The  state  estimates  generated  by  these  Kalman 


I filters  are  then  combined  using  a weighted  sum  with  the  a pos- 

teriori hypothesis  probabilities  as  weighting  factors.  If  one 
of  the  selected  parameter  vectors  coincides  with  the  true  para- 

r 

meter  vector,  this  algorithm  gives  the  minimxim  variance  state 
and  parameter  estimates.  Algorithms  for  filtering,  smoothing, 
and  prediction  are  derived  for  linear  and  nonlinear  systems. 

I They  are  described  in  a tutorial  fashion  with  results  stated 

explicitly  so  that  they  can  be  readily  used  for  computer  imple- 
I mentation.  Approaches  for  the  extension  of  MMEJ^  to  a more  general 

: class  of  adaptive  estimation  problems  are  outlined.  Several 
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1.  INTRODOCTIOM 


During  the  past  decade  considerable  advances  have  been 
made  In  the  theory,  algorithms,  and  applications  of  stochastic 
estimation  problems  Involving  linear  and  nonlinear  dynamics.  Th# 
linear  Kalman  filter  [1]  and  Its  diverse  extensions  to  the  nonlin- 
ear case  [2,3,4]  are  well  established  theoretical  and  algorithmic 
tools  with  extensive  applications. 

In  roost  practical  applications  of  recursive  estimation 
theory,  there  are  difficulties  In  obtaining  an  exact  mathematical 
model  of  the  physical  dynamic  process.  The  uncertain  parts  of  the 
syatetb  are  sometime  represented  by  an  unknown  parameter  vector. 
Examples  of  this  kind  include  the  ballistic  coefficient  and  lift- 
ing parameters  modelled  In  the  dynamics  of  a reentry  vehicle 
[4,5,6,7,81.  When  the  state  estimation  for  this  type  of  system 
has  to  be  carried  out,  the  variations  of  these  parameters  and 
their  Identification  play  a critical  role. 

Many  approaches  have  been  proposed  in  attempting  to 

* 

perform  state  estimation  together  with  parameter  identification. 
One  very  well-known  on-line  Identification  method  is  to  model  the 
unknown  parameter  as  a Markov  process  with  variance  related  to 

References  in  this  category  are  too  many  to  list,  one  may  consult 
the  IEEE  Transactions  on  Automatic  Control  (Dec.  1974) , a special 
Issue  on  system  Identification,  and  reference  [9]  for  listing 
of  related  references. 
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the  system  structure  and  the  range  of  parameter  variation.  The 
restriction  of  this  method  is  that  its  performance  is  critically 
Influenced  by  the  system  structure,  parameter  variation,  and  the 
required  bias  and  random  errors.  This  technique  usually  works 
well  within  a rather  small  region  of  the  state  apace  and  the 
variance  of  the  process  noise  can  only  be  determined  by  engineer- 
ing intuition  and  extensive  simulation  study.  This  method  how- 
ever, has  been  able  to  produce  excellent  estimation  accuracies 
in  reentry  vehicle  tracking  applications  [5,6,8]. 

There  exists  an  adaptive  filtering  and  parameter  identi- 
fication method,  which  we  shall  call  Multiple  Model  Estimation 
Algorithm  (MMEA)  in  this  report,  which  has  attracted  considerable 
attentions  in  the  academic  field  [10,  11,  12,  13,  14],  This  algor- 
ithm was  first  introduced  by  Magill  [10]  and  later  refined  by 
Lainlotis  [11] . The  estimation  algorithm  was  extended  to  adaptive 
control  by  Willner  [12]  and  Upadhyay  and  Lainlotis  [13] . 

The  basic  concept  of  NNEA  is  to  construct  a bank  of  Kalman 
filters  with  each  matched  to  a possible  parameter  vector  value. 

The  state  estimates  generated  by  these  Kalman  filters  are  then 
combined  using  a weighted  sum  with  the  posteriori  hypothesis  prob- 
abilities as  weighting  factors.  If  one  of  the  selected  parameter 
vectors  coincides  with  the  true  parameter  vector,  this  method  gives 
the  minimum  variance  estimates  of  both  the  state  vector  and  the 
parameter  vector-  In  most  physical  problems,  one  usually  has  a 

A 
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good  Idea  of  the  possible  values  that  a parameter  may  attain. 
Furthermore*  the  construction  of  the  MMEA  with  a steady  state  Kal- 
man filter  bank  requires  only  moderate  computation.  It  therefore 
has  attracted  some  attention  for  real-time  applications  tlS*  16]. 

The  purpose  of  this  report  is  to  Introduce  the  Multiple 
Model  Estimation  Algorithm.  It  will  be  described  in  a tutorial 
fashion  with  results  stated  explicitly  so  that  they  can  be  readily 
used  for  computer  implementation.  Furthermore*  the  discussions 
on  prediction  and  smoothing  are  believed  to  be  new.  Only  the 
algorithms  for  discrete  time  system  will  be  discussed.  This  is 
because  that  the  modern  estimation  and  control  algorithms  are 
mostly  implemented  on  digital  computers.  Due  to  the  fact  that 
MMEA  is  theoretically  more  sound  than  the  previous  methods,  it 
may  be  a potential  candidate  in  trajectory  re-construction  appli- 
cations. 

This  report  is  organized  as  follows.  In  the  next  section, 
the  problem  of  state  estimation  with  unknown  parameters  is  form- 
ulated. Possible  solutions  are  discussed  in  a tutorial  fashion. 

In  section  three,  the  Multiple  Model  Filtering  Algorithm  (MMFA) 
is  derived.  The  extensions  to  prediction  (MMPA)  and  smoothing 
(MMSA)  are  presented  in  section  four.  Discussions  of  the  first 
four  sections  assiune  linear  system  and  measurement  equations. 

The  extension  to  the  nonlinear  system  and  methods  of  algorithm 
realization  are  presented  in  section  five.  A simple  second  order 
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•xftnqpls  is  Includsd  In  ssction  six  to  illustrate  the  theory. 
Discussions  are  vfiven  in  the  last  section.  Two  appendices  which 
list  the  linear  smoothing  algorithms  and  the  Kalman  and  the  ex- 
tended Kalman  filter  algorithms  are  included  for  the  reference 
purpose. 
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2.  PROBLEM  FORMULATION 
2 . 1 Introduction 

Consider  a linear  stochastic  dynamic  system  whose  dynam- 
ics depend  on  a parameter  vector  Let  us  write  its  equations 
in  the  standard  state  apace  representation  and  in  the  discrete 


time  case. 


State  Dynamics 


x(t  + X)  « k(x)x{t)  + B(^)u(t)  + L(]r)C(t) 


(2.1) 


Measurement  Equation 


2(t)  “ C(Y)x(t)  + 0(t) 


(2.2) 


Next  we  define  the  different  variables  associated  with  eqR.  (2.1) 


and  (2.2). 


The  scalar  t is  a discrete  time  index 


t Of  If  2f  ... 


(2.3) 


The  state  vector  x(t)  c is  an  n-dimensional  vector.  The  input 
or  control  vector  u(t)  e is  an  m-dimensional  vector.  The 
plant  noise  vector  ^(t)  t Rp  is  an  p-dimensional  vector.  We 
assume  that  ^(t)  represents  a zero  mean  discrete  white  noise 
sequence  with  )cnown  covariance  matrix  H(t)  - pxp  matrix  - i.e. 


E { ^(t)  } = O for  all  t 


(2.4) 


cov  ( £(t);C(T)  ) = E { C(t)^'^(T)  } = H(t)6(t,T) 


(2.5) 


where  6(t,t)  is  the  Kroenecker  delta 


6(t,T)  = 


1 iF  t = T 


0 iF  t T 


(2.6) 


Note  that  the  plant  noise  covariance  matrix  H(t)  is  symmetric 
and  at  least  positive  semideninite 


H(t)  = H'^(t)  5^  0 


(2.7) 


The  measurement  noise  vector  £(t)eR^  is  an  r-dimensional  vector. 
We  assume  that  £(t)  represents  a zero  mean  discrete  white  noise 
sequence  with  known  covariance  matrix  e(t)  - an  rxr  matrix  - i.e. 


E { 6 (t)  } = 0 


(2.8) 


cov  ( 6(t);0(T)  ] = E { 0(t)0NT)  } = 0(t)6(t,T) 


(2.9) 


G(t)  = 0Mt)  > 0 


(2.10) 


Furthermore  we  assume  that  the  plant  driving  noise  £(t)  and  the 
measurement  noise  ^(t)  is  independent  for  all  values  of  t and  x, 
i.e.  , 


cov  f C(t);0(i)  ) = 0 for  all  t,x 


(2.11) 


The  above  fix  the  dimensions  of  the  different  matrices  that 


appear  in  eqs.  (2.1)  and  (2.2).  Thus 
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A(y;  is  an  nxn  matrix 

B(y)  is  an  nxm  matrix 

L(y)  is  an  nxp  matrix 

C(0)  is  an  rxn  matrix 

2.2  The  Parameter  Vector  y 

We  have  explicitly  shown  the  dependence  of  the  state 
dynamics  and/or  of  the  measurement  equation  upon  the  parameter 
vector  We  assume  that  the  parauneter  vector  xeRg  is  a q-dim- 
ensional  vector  whose  elements  represent  the  key  parauneters . 

The  elements  of  the  parauneter  vector  x general 

known  only  approximately.  The  degree  of  accuracy  by  which  the 
elements  of  x ®^e  known  are  strongly  dependent  upon  the  accur- 
acy of  modelling  a physical  process  by  eqs.  (2.1)  and  (2.2) 
and  the  experiments  that  have  been  carried  out. 

In  general,  before  the  initiation  of  any  real  time  es- 
timation and/or  control  experiments,  i.e., prior  to  time  t=0, 
one  has  some  idea  of  the  nominal  value  of  the  parameter  vector, 
denoted  by  x^'  degree  of  uncertainty  (e.g. , standard 

deviations)  associated  with  the  nominal  parameter  values. 

For  the  above  reasons,  it  is  reasonable  to  view  the 
parauneter  vector  as  a.  random  vector.  All  prior  information 
about  Y can  be  captured  in  its  prior  probability  density  function 
which  we  shall  denote  by  p(Y).  At  the  very  least,  our  best 
guess  about  x»  prior  to  any  additional  real  time  experimentation. 


is  the  nominal  value  ^ which  we  can  view  as  the  unconditional 
prior  mean 


E { X > = 


(2.12 


The  degree  to  which  we  "believe"  the  nominal  value  can  be 
communicated  to  the  mathematics  by  the  prior  covariance  matrix 
r - a qxq  matrix  - of  y,  i.e. 


( 1 ; 1 1 = E { (1  “ 1^)  (l  - } A 


(2.13 


It  is  also  reasonable  to  assume  that  the  uncertainty  associated 
with  the  parameter  vector  ^ has  nothing  to  do  with  all  other  un 
certainties.  Thus  we  make  the  assumption 


x(o) , ^(t) , and  £(t)  are  independent 
for  all  values  of  t and  t 


(2.1 


2.3  The  role  of  Y in  Filtering  Problems 


First  of  all  let  us  consider  the  filtering  problem  in 
the  context  of  state  estimation.  To  be  more  precise  let  us  de- 
note by  the  symbol  Z(t)  the  total  measurements  obtained  from  the 
initial  time  T=0  to  the  present  time  t.  These  measurements  in- 
clude both  the  Inputs  applied  to  the  system  and  the  actual  noisy 
sensor  measurements.  Thus  if  we  assume  that  the  first  sensor 
measurement  is  carried  out  at  t=l,  and  that  the  first  input  is 
applied  at  t“0,  then  the  data  set  Z(t)  is  defined  as  follows 
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2(t)  - { 2(1),  z(2),  2(t),  u(o),  u(l),  U(t-l)}(2.15) 

In  the  state  estimation  version  of  the  filtering  problem  one  is 
interested  in  obtaining  in  real-time  a "good"  estimate  of  the 
actual  value  of  the  true  state  vector  x(t)  based  upon  the  avail- 
able data  set  Z (t) ; this  state  estimate  is  commonly  denoted  by 

x(t/t)  (2.16) 

and  the  state  estimation  error  is  denoted  by 

x(t/t)  A x(t)  - x(t/t)  (2.17) 

We  can  now  have  several  cases,  depending  upon  the  relative  uncer- 
tainty associated  with  the  parameter  vector  x* 

Case  1 Parameter  vector  known  exactly 

This  is  an  unrealistic  case  and  corresponds  to  the  random  vector 
X having  zero  covariance 

r = 0 (2.18) 

— o — 

so  that 

X=Xo  <2.19) 

Under  these  assumption,  and  the  further  assumption  that  all  other 

random  vectors,  namely 

x(o) , ^(t) , e (t) 
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are  Gaueeian,  then  the  standard  discrete  time  Kalman  filter  [1] 
generates  the  optimal  estimate  of  the  state  in  the  sense  that  the 
state  estimate  x(t/t)  is  the  true  conditional  mean  of  the  state 

x(t/t)  - E { x(t)/Z(t)  } (2.20) 

In  addition  one  can  calculate  of f-line»  again  through  the  discrete 
time  Kalman  filter  algorithm  the  true  conditional  covarianca 
matrix  £(t/t) 

E(t/t)  = cov  t x(t)  ; x(t'/Z(t)  1 (2.21) 

Case  2 Parameter  Uncertainty  relatively  '* small” 

In  this  case,  we  assume  that  the  actual  value  of  the  parameter 
vector  X "very  close"  to  its  nominal  value.  Thus,  in  this  case, 
the  parameter  vector  covariance  matrix  is  small. 

I I I I = small  (2.22) 

An  alternate  way  of  characterizing  this  is  by 

I i II  « II  H(t)  II,  II  II  « II  0(t)  II  (2.23) 

which  means  that  the  parameter  uncertainty  is  much  smaller  than 
the  uncertainty  induced  in  the  state  by  the  plant  noise  ^(t),  and 
the  errors  introduced  in  the  sensors  by  the  measurement  noise  6^(t)  . 
Under  these  circumstances,  one  can  usually  trust  the  robustness 

The  discrete  Kalman  filter  algorithm  is  stated  in  the  Appendix  A. 
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of  the  Kalman  filter,  as  described  in  Case  1,  to  still  generate 
"good"  state  estimates  in  the  sense  that 

ft(t/t)  - E I x(t)/Z(t)  } (2.24) 

I(t/t)  5 cov  ( x(t)  I x(t)/Z(t)  1 (2.25) 

Case  3 Parameter  Uncertainty  Moderately  low 

i llol I increases,  the  errors  of  modelling  the  true 
values  of  the  parameter  vector  Y by  its  nominal  value  ^ become 
more  significant  and  the  performance  of  the  standard  Kalman 
filter  begins  to  deteriorate.  In  this  Intermediate  case,  and 
especially  when  the  major  effect  of  the  parameter  uncertainty 
are  reflected  in  the  state  dynamics  (2.1),  rather  than  the  mea- 
surement equation  (2.2),  there  have  been  several  cures  that  have 
been  suggested. 

The  basic  rationale  is  that  the  Increased  parameter  in- 
certainty in  the  system  dynamics  causes  errors  in  the  calculation 
of  the  one-step  predicted  estimate,  ft(t  + 1/t) , of  the  standard 
Kalman  filter  algorithm.  These  errors  can  only  be  corrected  by 
paying  more  attention  to  the  measurements,  which  although  noisy, 
still  contain  "good"  information  about  the  true  state.  Techni- 
cally, this  can  be  accomplished  by  increasing  the  magnitudes  of 
the  gains  of  the  Kalman  filter  and,  hence,  the  bandwidth  of  the 
Kalman  filter. 

One  way  of  accomplishing  this  objective  is  to  artificial- 
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Iv  incr»M<  Mlsot^d  of  tho  plant  noiao  oovariaaoo 

matrix  H(t).  Thla  trick  haa  baan  oftan  rafarrad  to  aa  introduo- 
ing  faka  ahlta  noiaa.  If  ona  can  gat  away  with  lt»  in  tha  aanaa 
that  tha  atata  aatlmatlon  arrora  x(t/t)  ramaln  aeoaptably  amall, 
than  thla  procadura  la  daalrabla  bacauaa  ona  can  atill  complata 
tha  (paaudo)  covariance  matrix  £ (t/t)  and  tha  Kalman  filter  galna 
off-line.  However,  thla  proceaa  of  turning  the  Kalman  filter  la 
more  of  an  art  than  a science. 

The  same  philosophy  of  changing  the  magnitude  of  the 
plant  noise  covariance  matrix  £ (t) , but  on  an  on-line  "adaptive" 
mode,  is  by  monitoring  the  behavior  of  the  residuals  of  the  Kalman 
filter.  The  residual  vector  of  time  t,  r (t/t) , is  defined  as  the 
difference  between  the  actual  measurement  at  time  t,  s(t),  and 
the  predicted  measurement 

r(t/t)  ^ z(t)  - C(x>x(t/t  - 1)  (2.26) 

In  the  case  of  no  parameter  uncertainty  (£^  *■  O)  the  residuals 
are  known  to  be  zero-mean  white  and  their  covariance  matrix,  de- 
noted by  S(t/t),  can  be  calculated  from  £(t/t) . As  the  parameter 
xincertainty  increases  this  is  reflected  in  the  nature  of  the  res- 
iduals, in  the  sense  that 

(a)  biases  can  be  observed  i.e., 

E I r(t/t)  } / 0 (2.27) 
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(b)  th«  residuals  bnooiM  oorr«lat«d  in  tlm«r  so  that 


thay  caaaa  to  ba  a urtiita  nolaa  taquanca. 

A variaty  of  inathodi  that  carry  out  raal  tina  taats  In  tha  real- 
duals  and  subaaquantly  changa  on-lina  tha  alananta  of  tha  plant 
nolaa  covarlanca  matrix  can  ba  suggaatad.  One  of  tha  almplaat  to 
InplMMint  la  tha  ona  suggaatad  by  Jazwlnskl  [2,17).  Tha  prica 
that  one  pays  In  these  adaptive  filtering  methods  is  increased 
real-time  computations  associated  with 

(a)  real-time  testa  and  confutations  involving  the 
residuals 

(b)  subsequent  transformation  of  the  residual-derived 
information  into  changes  in  the  covariance  matrix 
S(t) 

(c)  on-line  calculations  of  the  covariance  equation  and 
of  the  Kalman  filter  gain  matrix 

From  a pragmatic  point  of  view,  these  adaptive  filtering 
algorithms  change  in  a time-varying  way  the  gains  and  the  band- 
width of  the  Kalman  filter,  as  modelling  errors  become  significant 
and  diagnosed  in  the  residuals.  If  well  designed,  they  can  be 
effective  in  adjusting  the  bandwidth  of  the  Kalman  filter. 

It  should  be  noted  that  there  is  a tradeoff  associated 
with  hlgh-galn,  high-bandwidth  Kalman  filters.  Hlgh-galn  Kalman 
filters  tend  to  decrease  mean  errors  rapidly;  on  the  other  hand 


pow«r  to  past  through*  and  this  can  causa  increased  RMS  errors  in 
the  estimates.  The  successful  prior  timing  and/or  adaptive  filter- 
ing algorithms  have  to  take  explicitly  into  account  these  mean 
errors  vs.  RMS  errors  tradeoffs. 

Case  4 Moderate  Parameter  Uncertainty 

As  the  parameter  covariance  matrix  increases  further* 
the  off-line  or  on-line  turning  of  the  basic  Kalman  filter  cannot 
be  counted  upon  to  produce  good  estimation  accuracy.  This  is  due 
to  the  fact  that  the  contributions  of  the  parameter  errors  to 
model  uncertainty  can  no  longer  be  taken  care  of  as  equivalent 
white  noise. 

In  such  cases*  one  has  to  increase  the  real  time  com- 
plexity of  the  algorithm  so  as  to  explicitly  carry  out  some 
on-line  parameter  estimation.  In  other  words*  in  order  to  be 
able  to  obtain  reliable  state  estimates*  one  has  to  obtain  better 
estimates  of  the  parameter  vector  ^ based  upon  the  real  time 
measurements.  In  other  words,  the  filtering  algorithm  has  to 
simultaneously  generate 

(a)  a state  vector  estimate*  x(t/t) 

(b)  a parameter  vector  estimate*  ^{t/t) . 

Unfortunately*  even  in  the  simplest  case*  the  joint 

state  and  parameter  estimation  problem  constitutes  a nonlinear 
filtering  problem.  It  is  well  known,  118]  to  [22]*  that  the  true 
optimal  solution  to  a nonlinear  filtering  problem,  in  the  sense 
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iiniritlTH  tm  tm«  oonditional  M«a  of  tho  atoto 
* f *(t)/S(t)  ) rcqulroa  tho  on^ltn#  aolutton  of  a aot  of  non* 
PQgtlol  difforontial  oquatlona  at  oach  ond  ovorv  tlao  o 
awos^yaint  la  ■odo.  For  «l»o«t  all  problMS  of  practical  In- 
portanca#  tha  raal  tina  ooa^utational  raaooroaa  foroa  tha  da** 
aigoar  to  uaa  a avdx^tiaal  filtaring  algorithm. 

Tha  alnplaat  suboptlnal  filtaring  algorithm  la  tha  ao- 
callad  axtandad  Kalman  filtar.*  A allghtly  nora  complax  algori- 
thm la  tha  ao-oallad  aaoond  ordar  [4]  or  gauaalan  [2,23]  filter. 

Tha  tachnlqua  that  la  uaad  to  daalgn  tha  axtandad  Kalman 
flltar  la  that  of  atata  auga»ntation.  Thua,  In  addition  to 
ag.  (2.1)  which  definaa  tha  dynamic  atochaatlo  evolution  of  the 
"natural"  n state  varlablea  one  wrltea  another  aet  of  difference 
equations  of  tha  form 


X(t  + 1)  • i(t)  (2.28) 

In  case  that  It  Is  Known  that  the  parameter  vector  ^ Indeed 
a constant.  If  the  parameter  vector  x Known  to  change  slowly 
with  tla>e.  then  the  simplest  way  of  modelling  this  is  by  the 
stochastic  difference  equation 

X<t  + 1)  - x<t)  + pvt)  (2.29) 

The  extended  Kalman  filter  algorithm  is  stated  In  the  Appendix  A. 
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^ttro  £(t)  !•  a "fakv**  t«ro  iMan  ^it«  nolsa  prooaas  with  oovar- 
ianoa  matrix 

cov  ( i»(t)  I ii(t)  1 ■ M(t)6(t»T)  (2.30) 

The  oovarianoa  matrix  N(t)  haa  to  be  suitably  aelected  by  the 
deaigner  to  reflect  how  rapidly  and  by  how  much  one  can  reaaon- 
ably  expect  the  parameter  x to  change  or  drift  from  ita  prior 
nominal  value.  We  remark  that  more  complex  dynamic  modela  than 
that  ahown  in  eq.  (2.29)  can  be  used  If  prior  information  on  the 
"dynamica”  of  the  parameter  x available.  The  extended  Kalman 
filter  algorithm  that  generates  the  state  estimate  x(t/t)  and  the 

A 

paraieeter  estimate  x(t/t)  haa  much  more  severe  computational  re*> 
quirements  than  the  algorithms  discussed  in  Case  3.  These  addl« 
tional  requirements  are  due  to  the  fact  that  at  each  measurement 
time  one  has  to 

(a)  update  an  (n  q)  - dimensional  vector*  the  number  (n) 
of  state  variables  plus  the  number  (q)  of  the  para- 
meters 

(b)  propagate  an  (n  q)x(n  q)  (pseudo)  covariance 
matrix  using  the  standard  extended  Kalman  filter  co- 
variance  propagation  formula. 

(c)  calculate  a new  (n  + q)xr  Kalman  gain  matrix 

We  remark  that  all  the  "tricks"  discussed  in  Case  3 which  involve 
the  prior  turning*  or  adaptive  turning  based  on  the  residual  be- 
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havlour»  can  ba  usad  In  this  case  also  to  change  the  "fake  white 
noise"  covariance  matrices  \(t)  and  M(t) . 

2.4  Discussion 

The  above  brief  semiphilosophical  discussion  points  up 
some  of  the  issues  associated  with  the  effects  of  uncertain  para- 
meters upon  estimation  problems.  One  can  visualize  the  "robust- 
ness" of  the  varying  complexity  Kalman  filters  described  in  Cases 
1 to  4 as  shown  in  Figure  2.1 

The  way  Figure  2.1  is  to  be  interpreted  is  that  if  the 
true  parameter  is  in  band  3,  then  the  estimatorf  discussed  in 
Cases  1,2  will  not  give  satisfactory  performance,  while  the  es- 
timators discussed  in  Case  3 will  give  good  estimates.  Needless 
to  say  the  relative  sizes  or  shapes  of  these  robustness  bands  are 
next  to  impossible  to  calculate. 

The  point  that  we  wish  to  stress,  is  that  if  the  true 
parameter  is  outside  the  robustness  band  4,  then  the  extended 
Kalman  filter  discussed  in  Case  4 cannot  be  trusted  to  generate 
good  state  estimates,  even  though  on-line  parameter  estimation  is 
accomplished.  The  basic  reason  for  this  is  that  the  covariance 
linearizations  associated  with  the  extended  Kalman  filter  become 
invalid. 

For  this  reason  we  shall  explain  in  the  next  section  how 
one  can  attack  the  problem  of  large  parameter  uncertainty  through 
hypothesis  testing  and  subsequently  suggest  a suboptimal  procedure 
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that  oaa  b«  Mad  for  problama  with  larga  per«Mt«r  unc«rtalnty» 
aa  wall  aa  auddan  tranaitiona  of  tha  paramatara  (as  It  la  tha  caaa 
vith  aanauvarln?  raantry  vahiclaa) . 


I 
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3.  MULTIPLE  MODELS  FOR  HYPOTHESIS  TESTING  AND  STATE  ESTIMATION; 

FILTERING 

3.1  Introduction 

In  the  previous  section  we  have  outlined  the  different 
methods  that  can  be  employed  to  carry  out  state  estimation  when 
the  system  dynamics  contain  uncertain  parameters.  We  have  con- 
cluded that  as  the  parameter  vector  variance  increases  one  is 
forced  to  employ  nonlinear  filtering  algorithms,  e.g.,  the  ex- 
tended Kalman  filter,  which  simultaneously  estimate  the  para- 
meter vector  and  the  desired  state  variables.  We  have  also  re- 
marked that  even  these  sophisticated  algorithms  will  break  down 
as  the  parameter  uncertainty  increases. 

In  this  section  we  present  the  next  most  obvious  level 
of  complexity  to  take  into  account  the  effect  of  uncertain  para- 
meters. The  first  and  simplest  case  is  to  subdivide  the  parameter 
space  into  regions  and  see  what  happens  to  the  state  estimation 
algorithm  when  such  a discreti2ation  of  the  parameter  space  is 
carried  out. 

3. 2 Discretization  of  the  Parameter  Space 

As  we  have  remarked  in  Section  2.2,  the  parameter  vector  ( 

I 

X is  a q-dimensional  vector.  In  most  physical  problems,  one  has  j 

some  prior  idea  of  the  physical  ranges  of  the  elements  of  the  para- 
meter vector  j.  This  engineering  knowledge  can  be  translated  into 

a subset  Q of  R ; the  physical  significance  of  n :.s  that  it  re- 
Y q T 

presents  all  reasonable  values  that  the  parameter  vector  can  ! 
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attain. 


The  next  step  la  to  select  a finite  set  of  parameter 
values  denoted  by 


Yi»  y 2'  * • * ' 


(3.1) 


These  parameter  vectors  are  scattered  in  the  region  n . 


3. 3 Towards  the  MMFA;  Assumptions 

Let  us  suppose  that  the  parameter  vector  y,  which  appears 
in  the  state  space  description  of  the  stochastic  dynamic  system 
(2.1)  - (2.2)  does  indeed  coincide  with  one  of  the  defined 
above.  However,  prior  to  making  any  measurements  we  do  not  know 
the  "true  index"  i. 

Needless  to  say,  the  above  assumption  is  not  true  in  any 
real  life  situation,  in  the  sense  that  the  true  pareuneter  vector 
X will  be  "near,"  but  not  identical  to,  one  of  the  l^'s.  Once 
more,  we  shall  postpone  discussion  of  this  issue  for  the  time  being. 

Under  the  assumption  that  indeed  ^ coincides  with  one  of 
the  questions: 

1.  What  type  of  an  algorithm  can  be  used  in  order  to 
generate 

a.  the  true  conditional  mean  of  the  state,  and 

b.  the  true  conditional  covariance  matrix  of  the 
state 

given  a set  of  past  measurements.  We  remark  that 
this  constitutes  the  standard  estimation  or  filter- 
ing question. 
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2.  What  type  of  an  algorithm  can  be  used  to  identify 
the  true  parameter  given  a set  of  past  measure- 
ments. We  remark  that  this  constitutes  an  identi- 
fication question. 

One  may  argue  that  in  many  applications  one  may  not  be  interested 

in  the  identification  question,  but  only  in  the  state  estimation 

* 

problem.  Nonetheless,  it  turns  out  that  one  cannot  answer  the 
questions  independently,  but  one  roust  obtain  the  answer  to  both 
questions  simultaneously. 

We  shall  next  formulate  the  problem  in  a mathematically 
precise  way,  and  then  summarize  the  solution  algorithm. 

3.4  The  MMFA;  Formulation 

For  each  value  of  y_,  denoted  by  let  us  redefine  the 
matrices  in  section  2 as  follows 


; i = 1,  2,  . . . , N 

We  remark  that  the  matrices  A^,  B^,  can  be  time-varying; 

their  time  dependence  is  not  explicitly  shown. 


In  the  context  of  tracking  RV's,  if  one  tracks  a ballistic  RV, 
and  the  ballistic  coefficient  is  viewed  as  the  uncertain  para- 
meter, then  one  is  usually  interested  in  both  state  estimation 
for  good  tracking,  and  parameter  estimation  for  discrimination. 
A similar  situation  exists  for  maneuvering  re-entry  vehicles; 
in  the  MaRV  case  one  is  interested  in  estimating  parameters 
that  are  characteristic  of  the  magnitude  and  direction  of  the 
maneuver  accelerations. 
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Using  the  above  notation,  one  has  a class  of  N distinct 


linear  stochastic  dynamic  systems  described  by 

State  Dynamics 

x(t+l)  - Aj^x(t)+Bj^u(t)+L^C(t)  ; i*l,2,...,N  (3.3) 

Measurement  Equation 

z(t)  - C^x(t)+e(t)  ; 1=1,2,. ..,N  (3.4) 

The  characteristics  of  the  Gaussian  plant  noise  ^(t)  are 
still  given  by  eqs.  (2.4)  « (2.7),  while  the  characteristics  of 
the  Gaussian  measurement  noise  ^(t)  are  still  given  by  eqs. 

(2.8)  - (2.11). 

In  addition  to  the  plant  noise,  measurement  noise,  and 
initial  state  uncertainty,  we  must  specify  the  parameter  vector 
uncertainty.  Under  our  assumptions,  the  random  vector  ^ can 
attain  a set  of  discrete  values  X2'  •••'  In'  view  of  this, 
X is  a discrete  random  vector. 

We  can  model  this  fact  by  a set  of  hypotheses.  Let  H 
be  a scalar  random  variable  ( a hypothesis  variable)  and  let 


H 


1' 


H 


N 


(3.5) 


denote  a set  of  events. 

The  interpretation  that  we  attach  to  the  event 


H = is  1 » Ij 
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and  w«  can  think  of  this  phanomnon  as  that  ** nature”  has  select' 
ed  the  j-th  linear  system,  defined  by  eqs.  (3.3)  and  (3.4)  and  has 
placed  it  inside  a black  box. 


Before  we  obtain  any  data  from  the  system  in  the  black 
box,  we  have  to  have  some  idea  of  the  prior  probability  of  vdiich 
system  is  in  the  black  box,  or  equivalently,  the  probability  that 
1 ■ Xi  each  i. 

Let  P^(0)  denote  the  prior  probability  that  the  i-th 
syst^  is  in  the  "black  box."  Thus 


P^(0)  A Frob(H-H^)  » Prob(x“l£)  (3.6) 


with 


P.  (0)  >0,  2-  P.  (0)  «=  1.  (3.7) 

^ i-1  ^ 


Thus,  the  probability  density  function,  p(K),  of  the  random  vari- 
able H is 

N 

p(H)  - E (0)6(H-H,)  (3.8) 

i»l  ^ ^ 


where  6(«)  is  the  Dirac  delta  function. 

Remark!  The  numerical  values  of  the  prior  probabilities  Pi(0) 
reflect  to  the  mathematics  our  best  guess  on  which 
models  are  more  likely  to  be  in  the  black  box  prior 
to  their  generating  any  data.  If  initially,  i.e.,  at 
time  t>^0,  any  one  of  the  models  is  equally  likely, 
then  we  would  select  the  P^(0)  by 
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^ i3,10) 

«nd  Mk«  a Mt  of  aoiao  Masuramits 

mil),  m{2),  ...,  2(t)  (3.11) 

fron  tha  ayttap  in  tha  black  box.  Aa  wa  have  done  in  Section  2 
ve  let  S(t)  danota  tha  aat  of  all  past  measurements 

• < u(0),  u(l),  u(t-l),  »(1),  *(t)  > (3.12) 

Define  the  pr<^abillties 

(t)  4 Prob (H-H^/Z (t) ) 

(3.13) 

■ Prob(x*Xi/*<^> ^ 

to  be  tha  probability,  given  tha  measurement  sat  S(t),  that  tha 
i-th  hypothesis  (i.a.,  tha  i-th  model)  is  the  correct  one. 

Clearly 

P^(t)  > 0 (3.14) 

N 

E P. (t)  - 1 (3.15) 

i-1  ^ 
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aivttn  all  of  th«  above  information  and  notation,  wa  can  list  all 
the  information  that  we  would  like  to  obtain,  ae  well  on  the  re- 
quired algorithms  to  conqputa  the  variables  of  interest. 

1.  The  conditional  mean  of  the  state 


x(t/t)  ^ E f x(t)/2(t)  ) (3.16) 


2.  The  conditional  state  covariance  matrix 


I(t/t)  A cov  t x(t)  ; x(t)/Z(t)  1 (3.17) 


3.  The  dynamic  evolution  of  the  posterior  proba- 
bilities Pj(t);  ideally  we  would  like  a 
recursive*~relation , i.e.,  P>(t+1)  can  be  com- 
puted from  the  (t) . ^ 

Remark:  The  conditional  mean  and  the  covariance  can  be  computed 

once  p(x(t)/Z (t) ) , the  true  conditional  density  function 
of  the  state  of  the  system  in  the  "black  box"  has  been 
obtained . 

3.5  The  MMPA:  Derivations 


We  shall  obtain  recursive  relationships  of  the  general 
conditional  density  functions  at  time  t-fl  given  at  time  t. 

We  start  by  evaluating  the  conditional  probability  den- 
sity function 


p(x(t+l)/Z(t+D)  (3.18) 

Use  of  the  marginal  density  yields 

p(x(t+l)/Z(t+D)  « Jp(x(t+1),  H/Z(t+l))dH  (3.19) 
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\ ■A'  . 


p^<t+l),H/lU+X))'^(X(t+l)/H,I(t+l))p(H/a(t+U  (3.20) 

pxtbldblllty  d«noity  p(fl/t(t4-i))  auk  h%  Mxltfn 
uainf  th*  iwtAtioa  of  oq.  (3.13)  a* 


p(H/l(t+D)  - (t-fl)  fi  (H-H^) 


(3.21) 


Subotltuto  oq*.  (3.20)  and  (3.21)  into  aq.  (3.19) » and  Intagrata 


to  oibtain 


p(je(t+l)/l(t+D)  - i Pj^(t+l)p(x(t+l)/H^,B(t+l))  (3.22) 

Rmarkt  Wa  knoa  that  tha  conditional  dansitiaa 

can  ba  qanaratad  by  a bank  of  N Kalman  filtara  whore  aach 
filtar  is  "matchad"  to  a distinct  model  * i,a.»i-th 
hypothesis . 

It  is  is^rtant  to  realize  from  basic  Kalman  filtering  theory  that 
tha  following  relationship  is  true  for  each  conditional  probabil- 
ity density 

p (z  (t+1) /H^  rx  (t+1) ) p (X  (t+1) /H^ » * (t) ) 

p(x(t+l)/H^,B(ti-l))  - — = p(z(t+l)/H.,ft(tn  T3.23) 


and  that 


p(*(t+l)/H.,Z(t))  - /p(x(t+l)/H.,x(t))p(x(t)/Hj^,Z(t))^(t) 

^ (3.24) 


27 


UtMoric*  (Mm  om  asaunptions  all  dansitlaa  appaarf aqa  • 
(3.23)  and  (3.24)  &ra  Qauatlan,  and  hanoa  thay  dan  ba 
^araotarlaad  by  thalr  Man  and  oovarlanoa  Mtrix. 

Tha  baaio  Idaa  it  to  oonatruot  a ban)c  of  M Kaliaan  flit- 

arat  aaeh  Kalman  f 11  tar  la  daalgnad  uaing  tha  apaoiflc  paraMtar 

matrioaa  L^»  E»  and  (tha  initial  atata  covar- 

lanoa  matrix).  Bach  Kalman  filtar  in  tha  bank  is  drivan  by  tha 

saM  input  saquanca»  u(t),  applied  to  th«k'  aystan  in  tha  "black 

boXf"  and  by  tha  actual  measuramant  aaquanca,  z{t),  generated  by 

tha  aystam  in  the  "black  box." 

Let  x^(t/t)  denote  the  state  estimate  generated  by  the 

i-th  Kalman  filter.  Hore  precisely,  x^(t/t)  is  defined  by. 

X^(t/t)  x(t)/H^,  2(t)  I -/x(t)p(x(t)/H^,  Z(t))^(t)  (3.25) 

Let  £j^(t/t)  denote  the  conditional  covariance  matrix 
associated  with  the  i*th  Kalman  filter.  More  precisely 

I^(t/t)  4 cov  ( x(t) ;x(t)/H^,S(t)  ) 

-E{  (x(t)-x^(t/t>)  (x(t)-Xj^(t/t))’’/Hj^,Z(t)  > 

• /(x(t)-Xi(t/t»(x(t)-x^(t/t))^  • p(x(t)/H^,Z(t))^(t)  (3.26) 

Remark t All  the  E^(t/t),  i-l,2,...N  are  precoroputable . 

In  essence,  from  each  Kalman  filter  mean  Xj^(t/t)  and  covariance 
Btatrix  Z^(t/t)  we  can  construct  the  Gaussian  density  function 
p(x(t)/H^,Z(t)). 
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Th«  n«xt  problwR  i«  to  gonerato  an  ovarall  aatlnate  of 
tha  stata»  x(t/t) » according  to  eq.  (3.16)  of  the  system  In  the 
"black  box."  In  addition,  it  is  helpful  to  generate  the  true 
error  covariance  matrix,  £(t/t) , according  to  eq.  (3.17),  so  that 
ve  have  an  idea  of  how  accurate  the  estimate  x(t/t)  of  the  true 
system  state  x(t)  actually  is. 

We  demonstrate  below  how  the  overall  estimate  x(t/t) 
can  be  generated  onue 

a.  The  individual  Kalman  filter  estimates  x. (t/t) 

are  available,  and  ^ 

b.  The  true  conditional  probabilities  P. (t)  de- 
fined by  eq.  (3.13)  are  available.  ^ 

Fr<xn  eq.  (3.22)  we  have 

N 

p(x(t)/2(t))-  2 P.  (t)p(x(t)/H,  ,Z(t) ) (3.27) 

i-1  “1  “ i 

x(t/t)  ■ E{x(t)/Z(t)}  - y*x(t)p(x(t)/Z  (t)  )dx(t) 

N 

- 'Z  P.  (t)fx(t)p(x(t)/H.  ,Z(t))dx(t) 

i»l  ^ ~ ^ ~ 

N 

- Z P, (t)x. (t/t)  (3.28) 

i-1  ^ 

Thus,  the  overall  state  estimate  is  the  probabilistically  weighted 
average,  by  the  posterior  (hypotheses)  probabilities  P^(t),  of  the 
state  estimate  generated  by  each  one  of  the  N Kalman  filters. 

To  derive  the  true  conditional  covariance  matrix  £(t/t) 
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proo««d  M follows I 


E(t/t)  4 cov  [x(t) jx(t)/E(t)l 

- s{(x(t)-x(t/t)  (x(t)-x(t/t))Va(t)> 

- /(X(t)-x(t/t))  (x(t)-x(t/t))’^p(x(t)/Z(t))^(t) 

(t)-x(t/t))  (x(t)-x(t/t))’'. 

• p(x(t)/H^,Z(t))^(t)  (3.29) 


After  some  algebra  we  obtain 
N 

E(t/t)  - £p.  (t)  (Z.  (t/t)-Mx.  (t/t)-x(t/t))‘ 
i=l  ^ ^ 1 - 

• (x^(t/t)-x(t/t))^l  (3.30) 

Note  that  ^(t)  cannot  be  precomputed  because  it  contains  the  real 
time  estimates  x^(t/t)  generated  by  the  Kalman  filters  in  addition 
to  the  posterior  probabilities  P^(t)  which  as  we  shall  see  require 
real  time  measurements.  The  only  remaining  problem  is  to  calculate 
dynamic  evolution  of  the  porbabilities 

P^(t)  - Prob(H-H^/Z(t) ) 

- Prob(x“lj^/Z  (t) ) (3.31) 


We  will  relate  each  Pj^(t-H)  to  the  P^(t)  and  other  quan- 
tities that  can  be  found  from  Kalman  filters.  The  interesting 
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Mpttot  of  this  OAloulAtion  it  that  « truly  rtourtlvt  rtlttion- 
•hip  ot;:i  b«  obttintd  rtltting  quantititt  only  at  tuccaativo  mtat- 
UTtfMnt  ti«Mt»  t and  t-fl*  with  ralativaly  tnall  oonputational 


liMAirdtn. 


Towards  this  qoal  wa  prooaad  as  follows.  Considar  tha 


conditional  danaity 


p(H/2(t+D)  • (t+l)«(H-H.) 

i-1  * ^ 


0.32) 


Usa  of  Bayas  rula  yialds 


p(H/2(t+D)  - p(H/8(t+l)  ,2(t)) 

_ p(H,i.(fO)/Z(t)) 
p (2Tt-fl)/Z  (ti ) 


0.33) 


p(H/Z(t))  - 53  Pi <t)6(H-H.) 

i=l  ^ ^ 


(3.34) 


Note  that  according  to  our  notation  Z(t+1)  ■ {Z (t) ,2 (t+1) } 
Equations  (3.32)  to  (3.34)  yield 


P^(t+1) 


p(z(t+l)/Hj,Z(t)) 

p(z(t+l)/Z(t)) 


P^  (t) 


(3.35) 


The  density  p(z  (t-t-l)/H^,Z  (t) ) is  Gaussian  and  can  be 
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oaloulatftd  from  tho  l>th  Kalman  filtar 


p(«(t+l)/H^,2(t))  N(C^(t+l)x^(t+l/t)  (3.36) 


whara 


Sj^(t+1)  - C^(t+l)Z^(t+l/t)C^(t+l)+G(t+l) 


(3.37) 


Nota  that  tha  quantity  (t+Dx^^  (t+l/t)  Is  tha  predictad  maasura- 
roant  at  tinve  t4-I  generated  by  the  1-th  Kalman  filter. 

The  matrix  S(t4>l)  is  the  residual  covariance  matrix  asso* 
dated  with  the  1-th  Kalman  filter.  Mote  that  the  residual  co- 
variance  matrices  S^(t-t-l)  can  be  calculated  off-line  for  each 
Kalman  filter. 

It  remains  to  calculate  the  density  p(z  (t'M)/Z(t) ) in 
eq.  (3.35).  Use  of  the  marginal  density  leads  to 


p(*(t+l)/Z(t))  »/p(z(t+l),  H/Z(t))dH 

- /p(z(t+l)/H,2(t))p(H/Z(t) )dH 

M 

- /p{*(t+l)/H,Z(t))  5^  P.  (t)6(H-H.)dH 

j-1  3 3 


j 


N 

E (t)p(2(t+l)/H.,Z(t) 

i-1  J ~ 3 


(3.38) 


Remark:  Once  more  all  the  densities  p (z (t+1) /H . , Z (t) ) are  avail- 

able from  the  bank  of  Kalman  filters;  see  eqs.  (3.36) 


The  notation  N(a,A)  denotes  a Gaussian  density  with  mean  a and 
covariance  A. 
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O.S7).  Substituting  uq.  (3.38)  into  oq.  (3.35) 
yiolda  the  dosirod  result  that  the  dynamic  evolution  of 
the  probabilities  Pj^(t)  is  given  by 


p(s(t-^l)/H.  ,l(t)) 

P^(t+1)  - -,5 9^(t)  (3.35) 

£ P.  (t)p(i(t+l)/H^,Z(t)) 
j.l  J “ 3 


where  if  we  recall  that 


r ■ dim  s(t)  ■ number  of  measurements  (3.40) 


then 


with 


r 

■J 


1 

"I 


p(s(t+l)/H^,Z  (t))  - (2trl  [det  S^(t+l)l  • 

•exp|-y(*(t+l)-Cj (t+l)x^ (t+l/t) )^  Sj^(t+1) 
• (z(t+l)-C^ (t+l)x^(t+l/t))| 


(3.41) 


Sj(t+1)  - (t+l)Zj  (t+l/t)C^(t+l)-i-e(ti'l)  (3.42) 

The  relation  (3.39)  becomes  more  transparent  if  we  introduce  a 
some%diat  sing)ler  notation. 

Let  us  define  the  residual  (an  r-dlmenslonal  vector) 
vector  generated  by  each  Kalman  filter  by 


r^(t+l)  ^ 2(t+l)-C^(t+l)x^(t+l/t;  (3.43) 
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i.e.,  the  difference  between  the  actual  measurement  and  the  pre- 
dicted measurement. 

Then  from  each  Kalman  filter  we  can  obtain  the  scalar 
quantity  in  real  time 

Wj^(t+1)  A r^(t+l)'^sT^(t+l)rj^(t+l)  (3.44) 

Also,  let  B^(t+1)  denote  the  scalar  precompu tab le  quantity 

_r  _1 

B^(t+1)  A [27rj  ^[det  S^(t+l)j  ^ (3.45) 

Using  the  above  notation,  the  conditional  density  (3.41)  can  be 
written  as 

p(z(t+l)/Hj,Z(t))  = 6^(t+l)  exp|-|wj^(t+l)J  (3.46) 

From  eqs.  (3.46)  and  (3.39)  we  can  now  write  the  dynamic  evolu- 
tion of  the  probability  density  function  as 

6i(t+l)exp|-|  W^(t+l)| 

P^(t+1)  - ^ P^(t)  (3.47) 

Wj(t+l)}p.  (t) 

The  above  formula  illustrates  that  all  measurements  up  to  time  t, 
Z(t),  aie  captured  in  the  posterior  probabilities 

P^(t),  ?2(t)  , ...,  P^(t)  (3.48) 

The  new  measurement  at  time  t+1,  z (t+1) , influence  all 

''4 


I 


N residual  vectors  associated  with  the  bank  of  Kalman  filters 
according  to  eq.  (3.43)  and  generate  scalars  Wj^(t+1),  i«l,2,.,.,N. 
This  then  can  be  used  to  update  the  probabilities 

Pj^(t+1),  p2(t-H),  Pj^(t+1)  (3.49) 

according  to  eq.  (3.47).  Thus,  this  represents  a true  recursive 
solution  to  the  problem  of  probability  updates. 

A block  diagram  illustrating  the  MMPA  is  shown  in  fig- 
ure 3.1. 

3.6  The  MMFA;  Parameter  Identification 

In  the  previous  subsection,  we  have  described  the  basic 
idea  of  the  Multiple  Model  Filtering  Algorithm.  In  addition,  we 
have  derived  algorithms  for  MMFA  realization.  In  this  subsection, 
we  will  show  that  the  MMFA  for  parameter  identification  can  be 
obtained  in  a straightforward  manner.  The  minimum  variance 
estimate  of  the  unknown  parameter  is  the  conditional  mean  i.e., 

X(t)  - Jx  p(l/Z(t))d]r  = E{y/ii(t)}  (3.50) 

Recalling  the  fact  that  the  events  and  equivalent, 

we  can  rewrite  eqn.  (3.21)  as 

N 

p(l/Z(t))  = ' (l-Xi)  (3.51) 

v/hore  P^(t)  is  interpreted  as  the  pi^obt  oility  that  true 


3-~ 


measurements 


multiple  model  filtering  algorithm  block  diagram 


based  upon  all  the  data,  Z(t).  Using  (3.51)  in  (3.50)  yields 

£(t)  - ^i^t)  (3.52) 

The  covariance  of  £ can  be  ootained  similarly.  Assuming  that  £ 
is  unbiased,  then 


* cov(£(t)) 


./<x  - - £(t))Tp(^/z(t))di 

» (t)  “ £(t))T) 


(3.53) 


3.7  Diacusslon 

We  now  discuss  the  asymptotic  properties  of  this  algor- 
ithm from  a heuristic  point  of  view.  If  the  system  is  subject 
to  some  sort  of  persistent  excitation,  then  one  would  expect  that 
the  residuals  of  the  Kalman  filter  associated  with  the  correct 
model,  say  the  1-th  one  will  be  "small**,  while  the  residuals  of 
the  mismatched  filters  (j?^i,  j**l,  2,  N)  will  be  "large". 

Thus,  if  i indexes  the  correct  model  we  would  expect 


W^(t)  <<  Wj(t)  for  all  j i (3.54) 

If  such  a condition  persists  over  several  measurements  equation 
(3.47)  shows  that  the  "correct"  probability  P^(t)  will  increase 
while  the  "mismatched  model"  probabilities  will  decrease.  To 
see  this  one  can  rewrite  (3.47)  as  follows. 
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(t+1) 


P^(t)  ■ J (t+l)«xp|-jW^  (t+l)|  ’ 

?1  (t)  [d-P^  (t) ) B^(t+l)#3qp|-^^  (t+X)} 

- L Pj(t)Bj(t+l)exp{^^(t+l)}] 


(3.55) 


Under  our  assun^tions 


exp  i(t)|  • ^ 
exp|-^^(t)^  • 0. 

Hence  the  correct  probability  will  grow  according  to 


P^(t+1)  - P^(t)  « 


P^(t) tl-Pi(t)lBi<t+l) 


Pj(t)8j(t+l)exp|-|Wj(t+l)| 


> 0 (3.56) 


which  demonstrates  that  as  P^(t)  -*>  1,  the  rate  of  growth  slows 
down. 


On  the  other  hand,  for  the  incorrect  model,  indexed  by 
j^i,  the  same  assumptions  yield 


-P. (t)P, (t)B. (t+1) 

P.  (t+1)  - P.  (t)  s J i * < 0 (3.57) 

Pj(t)6^(t+l)exp|-|w^(t+l)} 


so  that  the  probabilities  decreased. 
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Thtt  •' 


oonoluslons  hold  if  wo  rowrito  (3.47)  in  tho 


form 


?^(t+l)  - P^(t)  Pj(t)Bj(t*H)oxp|-J»l^(t+l)|j  ^ 

P^(t)  ^B^(t+l)oxp^-^^(t+l)| 

•Bj  (t+l)axp|-^^ 


(3.58) 


Tho  obovo  dlocuotion  points  out  that  this  "idontlflcatlon"  schema 
is  crucially  dapandant  upon  the  regularity  of  the  residual  behav- 
ior between  the  "matched"  and  "mis-matched"  Kalman  filters. 

As  pointed  out  in  reference  [16],  the  dynamic  evolution 
of  the  residuals  may  not  follow  the  above  regularity  assumptions. 
This  may  be  caused  by  errors  in  the  selection  of  the  noise  sta- 
tistics or  using  a steady  state  Kalman  filter  design,  among  oth- 
ers. To  be  specific,  suppose  that  for  a prolonged  sequence  of 
measurements  the  Kalman  filter  residuals  turn  out  to  be  such  that 


Wj^(t)  5 W2(t)  = ...  : Wj,(t) 


(3.59) 


Then 


exp|-^j(t)|  = a for  all  i 
Under  this  condition  and  using  (3.58),  we  can  see  that 
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that  ona  of  tha  aay  tha  is  dominant,  i.a.. 


$1^  > 8^  for  all  IfQc 


Xn  this  easa,  tha  right-hand  sida  of  (3.60)  will  ba  nagativa  for 
all  ifOc,  «dki^  BMsns  that  all  tha  P^(t)  will  dacraasa  whila  tha 
pr<^ability  P|^  (associatad  with  tha  dominant  Bj^)  will  incraasa. 
Tha  significanca  of  this  affact  is  that  tha  B^'s  ara  indapandant 
of  tha  rasiduals  and  thair  magnitudas  are  not  determinad  by  which 
modal  is  true.  This  issue,  which  has  not  bean  discussed  in  tha 
literature,  is  baliavad  to  tie  with  tha  "Idantlf lability"  ques- 
tion of  this  schema. 

Above  discussions  merely  point  out  possible  shortcomings 
of  this  sesame.  These  issues  may  ba  adequately  answered  if  wa 
could  address  tha  following  questions. 

(1)  A rigorous  proof  to  show  tha  asyng>totic  properties 
of  tha  hypothesis  probabilities.  To  tha  bast  of  our  knowledge, 
such  a proof  is  not  available  in  the  literature. 

(2)  How  would  the  hypothesis  probabilities  behave  if 
none  of  tha  models  coincide  with  the  true  model?  Moor  and  HawXas 
(14]  used  a distance  measure  to  show  that  tha  probability  associa 
tad  with  tha  modal  idiich  is  tha  closest  one  to  the  true  modal 
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will  oonvwr^  to  unity,  it  this  olaia  is  wsrrsntsd,  ons  nay  bs 
abis  to  dssign  an  ad^tiws  parsastar  discrstisation  schaaa  which 
ra-disoratisas  the  paraaatar  vector  within  the  paraaatar  aubspaca 
which  is  the  closest  to  the  true  aodal  as  datarainad  by  the  hypo- 
thesis probability  and  the  distance  aaasura. 

(3)  Answers  to  the  above  questions  will  certainly  shade 
light  to  the  identifiability  problm. 

Finally t let  us  rc*eaphasise  the  significance  of  this 
scheme  from  the  estimation's  point  of  view.  This  algorithm  is 
optimum  in  the  minimum  variance  sense  in  state  and  parameter  es- 
timation ^ the  discretised  parameter  space  indeed  contains  the 
true  parameter.  This  is  true  because s (1)  We  use  the  condition- 
al swan  as  the  estliaate  and  (2)  the  algorithm  was  derived  without 
using  any  approximations. 
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4,  hypotbssis  twthm  amp  stact 

SMOOWmiO  AMD  WttDICTIOM 

4.1 

Zn  th«  previous  Motion,  we  have  derived  the  multiple 
model  filtering  algorithm  for  state  estimation  when  the  systwa 
dynamics  contain  uncertain  parameters.  The  parameter  vector  is 
diseretised  to  cover  a range  of  physical  values  that  it  may  pos- 
sibly attain.  A bank  of  Kalman  filters  Is  built  with  each  match- 
ing to  a parameter  vector.  The  a posteriori  probability  of  a 
given  model  being  true  is  used  to  combine  the  output  of  these 
filters.  Algorithms  for  state  estimation  and  parameter  identi- 
fication are  derived. 

In  this  section,  the  multiple  modes  smoothing  and  pre- 
diction algorithms  (MMSA  and  MMPA)  are  derived. 

4.2  The  HWSA  and  MMPA:  Assumptions 

The  system  equations,  measurement  equations,  parameter 
space,  and  hypothesis  probability  assumptions  made  in  the  section 
3.4  are  the  same  for  the  MMSA/MMPA  derivation.  We  only  modify 
the  variables  of  Interest  to  as  follows: 

1.  The  conditional  mean  of  the  state 

X(T/t)  A E { x(T)/Z(t)  } (4.1) 

2 . The  conditional  state  covariance  matrix 

£(x/t)  ^ cov  [ X (t) ;x  (t)/Z  (t)  J (4.2) 


3.  Th«  dynwde  «volution  of  tho  posttrior  prob«billtlo« 

F^(t)r  09«io«  w«  would  liko  a rocursivo  rolation. 

Ilpiliirtttit  <1)  idion  T>tf  it  is  osllsd  prsdiotion. 

sbsA  T<t»  it  is  osllsd  ssnothing. 
whsn  T«tt  it  is  osllsd  filtsring  snd  this  part  of 
algorithm  has  alrsady  bssn  prsssntsd. 

<3)  Ths  conditional  msan  and  the  covariance  can  be  com" 
putsd  once  the  conditional  density  function  has  been 
specified. 

In  ths  following,  ws  re-state  various  forms  of  prediction  and 
smoothing  in  terms  of  the  evolution  of  p(x (t)/Z (t) ) . 

(1)  Fix  lag  prediction/smoothing:  update  p(x(T)/S(t)) 

from  p(x(T-l)/2 (t-1) ) where  T-t  is  a fixed  constant 

(2)  Fix  interval  prediction/smoothing:  update 

p(x(x)/Z(t))  from  p(x(T-l)/Z  (t) ) 

(3)  Fix  point  smoothing:  update  p(x(T)/Z(t))  from 

p(X(T)/Z(t-l)) 

4 . 3 The  MMSA  and  MMPA;  Derivations 


Similarly,  we  start  by  evaluating  the  conditional  prob- 
ability density  function 


p(x{T)/Z(t) ) 

Using  the  marginal  density  yields 

p(x(T)/Z(t))  -^(x(t)  ,H/2(t))dH 

- ^(x(T)/H,Z(t))p(H/Z(t))dH 


(4.3) 


(4.4) 
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HoMV«r 


and 


p(H/Z(t)) 


t 


(3.21) 


P^(t)  - Prob (H-H^/2 (t) ) (3.13) 

Notice  that  P^(t)  la  Interpreted  as  the  probability  of  the  event, 
being  true  conditioned  upon  all  the  measurements,  Z(t). 

Unlilce  the  state  and  the  covariance  ((4.1)  and  (4.2)).  The  hy- 
pothesis probability  is  only  a function  of  one  time  variable,  i.e. , 
the  time  index  of  the  measurement  space.  Using  (3.21)  and  (3.13) 
in  (4.4)  yields 

N 

p(x(T)/Z(t))  • E P.  (t)p(x(x)/H.  ,2(t))  (4.5) 

i-1  i i 

This  equation  is  analogous  to  equation  (3.22).  Using  (4.5),  we 
obtain  the  predicted/smoothed  state  and  covariance  as 

x(r/t)  - E { x(T)/Z(t)  } 

■yx(T)p(x(x)/Z  (t)  )dx(x) 

N 

■ S P, (t)x. (x/t)  (4.6) 

i»l  ^ 


Zir/t)  » cov  (x(x)  ; x(x)/Z(t)l 

N y.  _ 

= P4(t)/(x/x)  - x(x/t))(x(x)  - x(x/t))’^ 

i-1  ^ “ 
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^p(x(T)/Hj^,S(t))dx(T> 


■ P^it)  + (X^(T/t)  - X(T/t) 


• (Xj^(f/t)  - x(T/t))^l 


(4.7) 


whftr«  Xj^(T/t)  la  the  estimate  from  the  i-th  smoother/predictor  and 
£^(x/t)  la  the  covariance  of  x^(T/t). 


Remarha:  (1)  The  realisation  of  NN8JV/NMPA  again  requires  a bank 
of  smoother/predictor  with  each  matching  to  a possible 
parameter  vector.  The  algorithms  for  the  individual 
smoother/predictor  realisation  have  long  been  made  avail- 
able,  for  example,  see  [3,  24-28],  or  Appendix  B. 

(2)  From  the  above  derivation,  the  hypothesis  probabil- 
ities (t)  for  «mx>thing/predietion~are  ^e  SMse  as 
those  for  filtering.  ¥he  dynamic  evaluation  of  t/Tt)  is 
•till  computed  by  using  equation  (3.47).  Recalling  that 
Pj^(t)  is  recursively  uj^ated  by  using  the  filter  resi- 
duals. Since  the  filtering  results  at  time  t are  ob- 
tained prior  to  any  prediction  and  smoothing  based  upon 
Z(t),  the  probabilities  Pj^(t),  1*1,  ....,  N are  always 
available. 


(3)  Prom  equations  (3.52)  and  (3.53),  the  parameter 
estimate  is  obtained  as  the  weighted  average  of  dlscre- 
tised  parameter  vectors.  Again,  there  is  only  one  time 
index  which  is  the  index  of  the  measurement  space.  The 
smoothlng/prediction  algorithm  for  the  parameter  esti- 
mate is  therefore  the  same  as  the  filtering  algorithm. 


In  sunnary,  we  state  the  following  procedure  for  apply- 
ing MMSA/MMPA. 

(1)  Compute  filtering  results,  i.e.,  obtain  x^(t/t), 
^^(t/t),  Pj^(t),  x(t/t) , and  £(t/t)  from  the  algorithms  of  the 
previous  section. 

(2. a)  For  prediction,  apply  the  individual  predictor  to 
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obtain  x^(t<fk/t)  and  ^^(t<t-k/t) , i.a.,  Itarata 

x^(t+l/t)  - A^x(t/t)  + Bj^u(t) 

and 

I 

k times  with  x^(t/t)  and  ^^^(t/t)  as  initial  conditions  where  k 
defines  the  discrete  prediction  time.  The  combined  estimate 
x(t+k/t)  and  covariance  ^(t+k/t)  are  obtained  by  using  (4.6)  and 
(4.7)  with  the  hypothesis  probabilities  (t)  the  same  as  those 
obtained  in  step  (1)  (filtering) . 

(2.b)  For  smoothing,  apply  the  individual  smoother  (see 
references  [24*281  or  Appendix  B)  to  obtain  x^(t-k/t)  and  Zj^(t-kA). 
The  combined  estimate  x(t-k/t)  and  covariance  E^(t-k/t)  are  obtain- 
ed by  using  (4.6)  and  (4.7)  while  the  hypothesis  probabilities 

(t)  are  constant  for  all  k and  equal  to  those  obtained  in  step  (jL) . 
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5.  MULTIPLE  MODEL  ESTIMATIQM  M«QORITHM  FOR  NONLIHEAR  SYSTEMS 

In  this  section,  the  MMEA  for  nonlinear  systems  is 
discussed.  From  the  previous  section,  It  Is  known  that  the 
smoothing  and  prediction  are  rather  straightforward  extensions  « 
of  filtering,  only  the  filtering  algorithm  will  be  emphasized 
here. 

Similar  to  the  linear  case,  we  define  the  following  non- 
linear system  and  measurement  equations  corresponding  to  the  i-th 
discretized  parameter  vector, 

State  Dynamics 

x(t+l)  • f(x(t),  u(t),  i(t),  x^)  (5.1) 

Measurement  Equation 

z(t)  • h(x(t)  , 6 (t)  , Xi^  (5.2) 

The  plant  noise  ^(t)  is  defined  by  equations  (2.4)  - (2.7)  and 
the  measurement  noise  6(t)  is  defined  by  equations  (2.8)  - (2.11). 

• j 

The  same  as  in  the  linear  case,  there  are  three  separate 
steps  in  the  multiple  model  estimation  procedure,  namely,  the  gen- 
eration of  Individual  state  estimates  matching  to  a given  para- 
meter vector  , the  evolution  of  the  hypothesis  probability  and  the 
combination  of  the  individual  estimates.  Let  each  steps  be  dis- 


cussed individually  below. 

(1)  It  is  well-known  that  the  realization  of  the  optimum 


state  estimation  for  systems  modelled  by  (5.1)  and  (5.2)  involves 
solving  a set  of  countably  Infinite  differential  equations  (18  - 
22).  It  is  therefore  practically  impossible  to  obtain  these  in- 
dividual optimum  estimates.  Suboptimum  filters  will  have  to  be 
used  to  construct  the  filter  ban)c,  i.e.,  to  produce  Xj^(t/t) 
approximately . 

(2)  The  equation  for  updating  the  hypothesis  probabil- 
ity is  stated  in  equation  (3.39) 


p(z(t+l)/Ht ,Z(t) ) 


£ P.  (t)p(2(t+l)/H.,Z(t) 

1»1  ^ J 


3.39 


In  arriving  at  this  equation,  no  assumption  was  made  on  which  type 
(linear  or  nonlinear)  of  systems  are  being  considered.  It  is 

I 

therefore  still  valid  for  nonlinear  estimation.  It  however,  can- 
not be  calculated  exactly  due  to  the  fact  that  the  exact  realiza- 
tion of  the  individual  a posterior  density  p(z (t+l)/H^,Z (t) ) can 

not  ,be  obtained,  it  can  only  be  evaluated  approximated  with  a 

* 

sub-optimal  filter  (such  as  the  extended  Kalman  filter  ) for 

computing  Xj^(t/t)  and  ^^(t/t). 

(3)  Assuming  that  the  optimum  individual  estimate 

£^(t/t)  and  its  covariance  ^^(t)  are  available,  the  optimum  state 
estimate  and  its  covariance  can  be  computed  by 

The  extended  Kalman  filter  equations  are  listed  in  the  Appendix  A. 


(3.28) 


X (t/t) 


N 

p^(t)Xi(t/t) 


Ut/t)  - E l£i(t/t)  + (X.  (t/t)  -x(t/t)) 

i»l  ^ -1  -1  - 

*(x^(t/t)  - x(t/t))^l  (3.30) 

Similarly,  in  order  to  realize  (3.28)  and  (3.30)  for  states  and 
measurements  represented  by  (5.1)  and  (5.2),  one  has  to  use 
suboptimum  filters  to  generate  the  individual  estimates  x^(t/t) 
and  (t/t) . 

Let  us  re-emphasize  that  equations  (3.28),  (3.29),  and 
(3.30)  are  exact  representations  for  the  solution  of  the  nonlinear 
estimation  problem  for  systems  modeled  as  (5.1)  and  (5.2).  In 
other  words,  the  a posterior  hypothesis  probabilities  evolution 
and  the  method  of  computing  the  combined  estimate  are  optimum  if 
each  individual  estimate  can  be  obtained  optimally. 

Numerous  s"l)optlmum  filters  have  been  proposed  for  non- 
linear estimation  (2,4,28-33].  The  most  popular  filters  are  the 
extended  Kalman  filter  and  the  second  order  filter  [2,4]  among 
others.  Especially,  the  extended  Kalman  filter  has  attracted 
considerable  attentions  for  practical  applications  (2-8].  The 
second  order  filter  can  generally  provide  improved  performance 
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than  the  ext  >«  K'  *in«n  filter  with  the  trade-offs  of  higher 
computational  jurden.  A comparison  of  various  nonlinear  filters 
may  be  found  in  [34,35].  All  these  filters  tiay  be  used  for  the 
MMEA  realisation.  A specific  selection  may  be  based  on  a partic- 
ular physical  problem  and  the  required  performance.  For  real-time 
application,  one  usually  favors  a simple  filter  pending  on  the 
available  con^uter  resources.  For  off-linear  processing  espe- 
cially in  the  post-mission  smoothing  application,  a sophisticat- 
ed algorithm  is  usually  preferred. 
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6.  WMPLE 

In  thin  section,  we  present  an  example  to  Illustrate 
the  theory.  Only  the  filtering  algorithm  is  tested. 

Consider  the  following  second  order  continuous  systMi. 


(6.1) 


This  system  may  be  used  to  describe  the  motion  of  a vehicle  along 
a given  axis  with  drag  (represented  by  ”y")  &nd  control  force 
(represented  by  "u")«  if  denotes  the  target  range  and  a radar 
is  used  to  take  range  measurements,  the  measurement  equation  is 


2 « Xj^  + n 


(6.2) 


where  n is  measurement  noise.  The  measurements  are  taken  at 
discrete  Instemce  of  times.  A corresponding  discrete  system  of 
(6.1)  Is 


. 

Xj(k+1) 

X2(k+1) 

^ m 

( 

Y 

e“YAt 


» e 

m ■ 

(k) 

0 

Xj(k) 

• « 

+ 

Y 

u 


(6.3) 


where  tt  Is  the  time  interval  between  measurements.  A multiple 
model  filter  Is  used  to  estimate  Xj^,  x^  and  to  identify  y.  Three 
Y values  are  assiimed,  i.e.,  y**0.,  .5,  or  1..  The  system  and  con** 
trol  matrices,  A and  L for  those  y values  with  sampling  Interval 
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tt  • ,1 


il)  y - 0. 


The  measurement  noise  standard  deviation  is  equal  to  one.  The 
time  initial  state  is 

- 100.  ; *2  " 

The  following  convention  is  used  to  relate  the  hypothesis  to  the 
parameter  values. 

— ► Y - 0. 

Kj  Y - .5 

Hj  — Y - 1. 


S2 


Two  •xporinonta  ar«  porfocMd.  Ihwy  «r«  doworibod  in- 
dividually balow. 

Bxpayiwant  It  Paranatar  y staya  constant,  control  d is  aqual  to 
taro. 

Thraa  casas  with  tha  trua  parainatar  baing  aqual  to  ona 
of  tha  thraa  poasibla  valuas  in  aach  casa  ara  tastad.  Tha  a 
postariori  hypothasis  probabilitias  for  all  thraa  casas  are 
plotted  in  Figure  6.1.  Tha  initial  hypothesis  prc^abilitias  ara 
uniformly  distributed.  Tha  trua  system  is  always  identified  in 
within  10  data  points 

Experiment  2:  Parameter  y jumps  between  models,  control  u is 

nonsero. 

Tha  control  force  is  assumed  to  be  equal  to  50  and 
known  to  the  estimator.  Assuming  the  Initial  time  is  zero,  the 
true  Y time  history  is 

Y * 0.  for  0 £ t .<  2 

Y • .5  for  2 < t < 4 

Y * 1.  for  4 < t < 6 

It  therefore  represents  a y history  with  sudden  jumps.  The  y 
estimates  are  presented  in  Figure  6.2.  Notice  that  the  filter  is 
always  able  to  identify  the  true  system.  Two  modifications  are 
implemented  in  the  algorithm  in  this  case. 

(1)  The  hypothesis  probabilities  are  hard  bounded. 

This  is  to  prevent  any  probabilities  from  converging  to  zero  (or 
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om)  . Hh«n  it  dots#  it  will  bm  vwry  diffioult  for  tho  probabili- 
tios  to  branch  out  again  whan  tha  trua  syatam  haa  actually  switch- 
ad*  Ifhrn  boui^  usad  in  this  axparisisnt  is  howavart  vary  saall#  i.a.« 

PrlB^/8(t))  > .0005  for  i-1,2,3 

(2)  Although  thara  is  no  process  noise  assumed  in  the 
system,  a process  noise  term  with  covariance 


is  used  In  the  filters.  This  is  included  also  for  the  purpose 
of  preventing  the  filter  from  being  over-confident  in  its  esti- 
mates therefore  not  able  to  switch  to  a different  system.  If 
there  is  no  process  noise  added,  the  estimates  of  a mis-matched 
filter  can  drift  far  away  from  the  true  states.  When  the  true 
parameter  jumps  to  a different  value,  i.e.,  an  originally  mis- 
matched filter  now  becomes  matched,  it  takes  extremely  long  per- 
iod of  time  for  the  algorithm  to  identify  the  true  system  again. 
Leaving  proper  process  noise  level  in  the  filter  will  keep  the 
mis-matched  filter  estimates  sufficiently  close  to  the  true  state 
so  that  the  algorithm  is  adaptive  to  the  parameter  jumps.  The 
control  variable  u also  plays  a critical  role  in  this  experiment. 
It  represents  a persistent  excitation  to  explore  differences 
among  these  systems.  A basic  issue  which  still  needs  answer  is 
on  the  input  design  for  system  identification  in  using  MMEA 
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The  above  two  experiments  are  simple  but  illustrative. 
The  first  experiment  indicates  that  the  MIFA  oan  quickly  identify 
the  true  system  with  a constant  parameter.  For  time-varying 
parameters,  some  modifications  are  necessary  so  that  the  algori- 
thm is  adaptive  to  sudden  parameter  changes. 
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7.1  SvMMry 

In  this  rnport,  w«  havn  discussed  the  problem  of  state 
eetlmatlon  with  uncertain  parameters  and  presented  the  solution 
by  utilisating  the  multiple  model  estimation  algorithm  (MNEA) . 

The  following  summaries  pertaining  to  the  properties  of  NMEA  are 
listed  without  any  specific  order. 

(1)  Theoretically,  the  MMEA  provides  the  minimum  vari- 
ance estimates  of  both  state  amd  parameter  if  one 
of  the  chosen  models  coincides  with  the  true  model. 

(2)  If  the  a posteriori  hypothesis  testing  probabilities 
converge  asymptotically,  the  true  parameter  is  iden- 
tified with  probability  one. 

(3)  The  hypothesis  probabilities  for  smoothing  and  pre- 
diction are  the  same  as  those  for  filtering. 

(4)  The  hypothesis  probability  update  equation  and  the 
weighted  sum  equations  are  optimum  in  the  minimum 
variance  sense  and  they  are  the  same  for  both  lin- 
ear and  nonlinear  systems. 

The  usefulness  of  MMEA  can  only  be  fully  understood  and 
evaluated  after  applications  to  significant  physical  problems. 
Applications  to  the  trajectory  estimation  area  have  still  to  be 
carried  out.  The  application  to  the  F-8C  airplane  real-time  con- 
trol system  (16]  has  shown  encouraging  results  and  suggested 
further  study  areas  in  theory  as  well  as  in  algorithm  design. 


7.2  Oiaouaciont  Bxtanaion  to  a Class  of  TlSM-Varying  Para* 
and  Suboptlaal  Amproaoiias 

Strictly  apaaXing,  tha  HHEK  prasantad  In  this  raport  is 
optisMB  only  for  systwM  with  tina* Invariant  paramatars.  Tha 
thaoratical  and  practical  iaplications  of  using  NMBA  to  systams 
with  tiiM*varying  paramatars  are  not  complataly  understood.  The 
exampla  in  tha  previous  section  has  clearly  indicated  that  some 
modifications  must  be  Incorporated  in  order  to  make  the  MNEA  to 
follow  parameter  jus^s.  This  is  because  once  the  true  parameter 
is  identified,  tha  algorithm  is  locked  on  tha  triae  system  and 
tha  mis*matched  Kalman  filter  begins  to  drift  away  from  the  true 
state.  When  the  true  parameter  has  switched  to  a different  value, 
it  usually  takes  a long  time  for  the  algorithm  to  branch  out 
again.  The  requirement  for  a time*varying  parameter  MMEA  is  to 
stake  the  mis*matched  output  still  sufficiently  close  to  the  true 
state  and  to  keep  the  hypothesis  probability  from  coming  too 
close  to  zero  (or  unity) . 

Thera  is  a trivial  extension  of  the  MMEA  to  a special 
class  of  time*varying  parameters.  Consider  the  parameter  space 
which  contains  N parameter  vectors  each  with  dimension  q,i.e., 

Rq  • ‘ 1 - li  ' i-l,.../N) 

At  the  time  t,  the  true  parameter  is  equal  to  At  the  next 

instance  of  time,  the  true  parameter  may  be  equal  to  any  param- 
eters in  Rg.  As  time  progresses,  the  true  parameter  is  changing 
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around  with  its  values  within  R^.  Defining  two  types  of  hypoth-* 

eses  by 

H.  (t)  " the  hypothesis  that  1 ■*  Xi  tine  * t is  true#  it 
^ is  therefore  a local  nypouesis 

H.  (t)  ■>  the  hypothesis  that  a giving  history  for  time  up  to 
^ t of  Y indexed  by  k is  true,  it  is  therefore  a global 

hypothesis. 

These  two  types  of  hypothesis  are  related  by  the  following  equa- 


tions 


Hj,(t)  - H.  (t)  • H.  (t-1)  • • H.  (1) 

*'t-l  *^1 

where  the  index  for  k^,  k^  is  1,  N,  the  index  for  k 

is  1,  and  9 denotes  the  ”and”  operator.  It  is  clear 

that  each  defines  a possible  sequence  of  x history.  With 

this  definition,  one  may  proceed  in  parallel  to  the  development 
of  this  report  to  obtain  a new  MMEA  for  time-varying  parameters. 
The  derivation  is  briefly  stated  below. 

1)  For  state  estimate  and  covariance 

* 

Let 

P^(t)  - Prob(H(t)  « H^(t)/Z(t))  (7.1) 

for  i*l,  ...,  N^.  It  is  trivial  to  show  that 
x(t/t)  »•  E(x(t)/Z(t)) 
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(7.2) 


- £ P^(t)x^(t/t) 
i-1 

£(t/t)  « cov  (x(t) , X(t)/I(t)l 

- £ P^(t)  (E^(t/t)  + (X^(t/t)  - x(t/t)l 

i-1 

•(x^(t/t)  - xCt/t))*^)  (7.3) 

where  x^(t/t)  - E(x(t)/H^(t) ,2 (t) ) 

Z^(t/t)  * cov  lx(t),x(t)  /H^(t),2(t)) 

2)  For  prohablllty  update 

Using  the  conditional  probability  relation  yields 

p(H(t+l)/Z(t+D)  » p(H(t+l)/H(t)  ,2(t+l))p(H(t)/Z(t>H))  (7.4) 

Using  Bay's  rule  on  p(H(t+l)/H(t) ,2(t+l) ) yields 


p(H(t+l)/H(t)  ,Z(t+D) 
p(2(t-^l)/H(t+l)  ,H(t)  ,Z(t) ) 
p(2(t+l)/H(t)  ,Z(t)) 


p(H(t+l)/H(t)  ,Z(t)) 

I' 


Define  the  following  probability  density  functionE; 


(7.5) 
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p(H(t+l)/H(t)  ,Z(t+l) ) - ^ P(H^(tfl)/H(t)  ,Z(t+l))6(H-HjL)  (7.6) 

N _ 

p(H(t+l)/H(t)  ,Z(t))  - P(H^(t+l)/H(t)  ,Z(t))6(H-H^)  (7.7) 

i=l 

where  P(H^ (t+l)/H(t) ,Z (t+1) ) = Prob(H(t+l)=H^ (t+l)/H(t) (t+1) ) 
and  P(Hj (t+l)/H(t) ,Z(t) ) is  the  probability  that  the  parameter 
will  switch  to  given  a past  history  of  and  all  the  past 
measurements.  It  is  determined  by  the  property  of  the  hypothesis 
process.  If  the  hypothesis  process  is  a Markov  process,  this 
probability  becomes  the  transition  probability,  i.e., 

P(H^(t+l)/H(t) , Z(t))  = P(H^(t+l)/H(t))  (7.8) 


= P(H^(t+l)/H(t)) 


For  exajnple,  if  the  parameter  may  change  to  any  parameter  in 
with  equal  probability,  we  may  assume 


P(H^(t+l)/H(t) , Z(t))  = I for  i=l,  ...,  N. 


Using  (7.6)  and  (7.7)  in  (7.5)  yields 
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P(H^(t+l)/H(t),Z(t+D) 

p(z(t-*-l)/H(t)  ,Z(t)) 

Using  the  equation 

p(if(t)/z(t))  - £ p^(t)«<S-H^)  <■'•10' 

k»l 


in  (7.9)  yields 

P(H^(t+l)/Hj^(t)  ,Z(t+D) 


p (z  (t+1)  (t+1)  (t)  pZ(t) ) 

p(z(t+l)/H,^(t)  ,Z(t)) 


P(H^(t+l)/H^(t) ,Z(t)) 

(7.11) 


where  p(2(t+l)/H^(t+.l)  ,H,^(t)  ,Z(t) ) is  the  residual  density  of 
the  filter  which  was  matched  to  the  Ic-th  history  and  is  now 

matched  to 

p(z(t+l)/H^(t) ,Z (t) ) 


* V p(z(t+l)/H_(t+l)  ,H.  (t)  ,Z(t))P(Hjjj(t+l)/H^(t)  ,Z(t)) 

(7.12) 

m=l 

Next,  we  relate  P (H. (t+1) /Z (t+1) ) to  the  conditional 
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probability.  Using 


P(Hj(t+l)/Z(t+l)) 

- y*  P(H^(t+l)/H(t)  ,2(t+l))p(H(t)/Z(t+l))dH(t)  (7.13) 

and  equation  (7.10)  one  obtains  the  following  equation  for  each 

Hj^(t) . 

P(H^(t+l),H^(t)/Z(t+D) 

» P (H^  (t+1)  /Hj^  ( t)  , Z (t+1) ) P (Hj^  (t)  /Z  (t+1)  ) 

for  i«l/  K and  )c«l,  (7.14) 

where  P(H^  (t+l)/Hj^(t)  ,2  (t+1) ) Is  specified  in  equation  (7.11). 
Notice  that  P (H^  (t+1)  ,Hj^  (t) /Z  (t+1) ) is  the  updated  hypothesis 
probability.  Next*  we  derive  the  equation  for  computing 
P(Hj^/Z  (t+1) ) . Using  Baye's  rule  on  P (Hj^/Z  (t+i) ) yields 

P(Hj^(t)/Z(t+l) ) 

p(2(t+l)/H.  (t)  ,Z(t)) 

- — P(H.  (t)/Z(t))  (7.15) 

p(z(t+l)/2(t)  ) ^ 

where  P(Hj^  (t)/Z  (t) ) is  the  a posteriori  hypothesis  probability 
at  time  ■ t,  i.e.,  Pj^(t).  The  probability  density  functions  of 
(7.15)  are  computed  by  using 
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p(s(t+l)/i^(t),S(t)) 

N 

■ £ P(»(t+l)/H,j(t),H^(t+l),E(t))P(H^(t+l)/gj^(t)  ,Z(t)) 
j-1  (7.16) 


And 


p(*(t+l)/Z(t)) 

- p(2(t+l)/Hj^(t)  ,2(t))Pjj^(t)  (7.17) 

Bt*l 

The  probability  update  is  therefore  carried  out  by  using  equations 
(7.11),  (7.14),  and  (7.15).  These  relations  can  be  further  con- 
densed with  the  following  simplified  notations. 

P(H^/Hj^,Z(t))  - P(H^(t+l)/Hj^(t)  ,Z(t))  (7.18) 

p(a(t+l)/H.  (t+1)  ,il  (t)  ,Z(t)) 

MHi/I,^)  - i 

5^  p(z(t+l)/H^(t+l)  »H^(t)  ,Z(t))P(Hj/Hj^,Z(t)) 
j-1  (7.19) 

a conditional  likelihood  ratio 
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pU(t+l)/^(t),l(t)) 

pCsCt-t-D/ECt)) 


MV  - 


M 

£ p(«(t+l)/V(t),H^(t+X),B(t))P(H^/V,«(t)) 



N 

E { E p(*(t+l)/f  (t),H^(t+l),E(t))P(H^/8_»*<t)))P«(t) 
m-1  n-1  ” « n n m 

(7.20) 

■ likelihood  ratio 

Using  (7.18),  (7.19),  and  (7.20),  equations  (7.11),  (7.14),  and 
(7. IS)  nay  be  combined  to  become 


Pj(t+1)  - l(H^/Hjj)P(H^/^,Z(t))MHj^)P^(t)  (7.21) 

Notice  that  £or  i«l,  ....,  N and  k«l,  ...,  M^,  the  index  for  j is 
1,  ....,  The  probability  update  is  carried  out  with  the 

conditional  probability  which  characterizes  the  hypothesis  process 
Itself  and  the  likelihood  ratios  which  use  the  new  information 
through  residual  density  functions  of  each  filter. 

The  MMEA  for  a constant  parameter,  i.e.,  the  algorithm 
discussed  in  Section  3,  is  only  a degenerate  case  of  (7.21). 

When  the  parameter  is  a constant,  the  local  hypotheses,  H^,  and 
the  global  hypotheses,  Kj^,  become  the  same.  The  number  of  hypothe- 
ses is  limited  to  the  nxunber  of  peurameter  vectors  in  R . Further- 

q 

more,  the  conditional  probability  of  equation  (7.11)  becomes 
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P(H^(t+l)/l^(t),Z(t)) 

• 1 wten  ®i  * * ®k 

■>  0 elae«rh«re 


Using  (7.22)  in  and  l(^)  yields 


when 


\ 0 elseidiere 

P(*(t+l)/H.  (t+1)  ,Z(t)) 

» = JS 

£ p(*(t-H)/H^(t+l),Z(t))P^(t) 


m«l 


(7.22) 


(7.23) 


(7.24) 


Using  (7.22) » (7.23),  and  (7.24)  in  (7.21),  we  obtain  equation 
(3.39),  the  probed>ility  update  equation  for  the  constant  parameter 
case.  This  completes  our  a posteriori  hypothesis  probability 
derivation . 

An  obvious  problem  with  this  algorithm  is  that  the  number 
of  Hj^(t),  is  growing  with  t.  In  order  to  malce  this  algo- 
rithm practical,  one  has  to  limit  the  growing  number  of  hypothe- 
ses. In  the  following,  several  subopt imal  approaches  for  the 
time-varying  parameter  MMEA  problem  are  outlined.  The  first  two 
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approach**  ara  aiaad  at  limiting  tha  uumbar  of  poaalbla  aaquanoaa 
(or  hypothaaaa) . Tha  laat  two  approachaa  ara  mainly  to  raduca 
tha  ^anca  of  tha  algorithm  baing  locked  on  a particular  aystam. 

(1)  Maximum  Likelihood  Probability  Approach 
Conaider  the  caae  that  at  time  t there  are  only  N 

hypotheaea  aelected.  For  the  next  time  period,  each 
hypotheaia  may  grow  with  M poaaibllitlea.  It  therefore 
haa  M • N hypotheaea  after  each  filter  update.  Theae 
M * N hypotheaea  are  then  limited  by  aelectlng  only 
thoae  M which  have  the  largeat  hypothesis  probabilities. 

(2)  Transition  Probability  and  Finite  Memory  Hypothesis 
Process  Approach 

Suppose  that  the  filtering  process  has  limited  mem- 
ory so  that  ^(t)  is  replaced  by  the  most  recent  local 
hypothesis  Hj^(t).  Furthermore,  it  is  assumed  that  the 
hypothesis  process  is  a Markov  process.  Then  one  is 
interested  in  updating 

P(H^(t+l)/Z(t+l) ) ; i-1,  N. 

from 

P(H^(t)/Z(t) ) for  all  k-1,  ...,  N. 

With  this  assumption  and  using  (7.13),  one  obtains 
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II 

P(E^(t+l)/l(t+D)  - £ P(H^<t+l)/H,^(t),S(t+l)) 

k-1 


• P(H^(t)/l(t+D)  (7.25) 

whttre  P(H^(t4>l)/Hj^(t)  ,Z(t<fl))  is  obtained  by  an  equa- 
tion similar  to  (7.11),  i.e., 

P(H^(t+l)/B,^(t)  ,a(t+D) 

p(s(t+l)/H.  (t-H)  ,H.  (t)  ,Z(t)) 

- — = i = P(H.  (t+D/H.  (t)  ,Z(t)) 

p(2(t+l)/Hjj(t)  ,Z(t))  ^ ^ (7.26) 

P(H^/Z(t-fl))  is  Obtained  by  an  equation  similar  to 
(7.15),  i.e., 

p(z(t+l)/H.  (t),Z(t)) 

P(IL/Z(t+1))  ■ 5 P(H.  (t)/Z(t))  (7.27) 

p(*(t+l)/Z(t))  ^ 


where 

p(z(t+l)/Hj^(t)  ,Z(t)) 

H 

» £ p(2(t+l)/Hj^(t)  ,H^(t+l),Z(t))P(Hj(t+l)/Hj^(t)  ,2(t)) 
j-1  (7.28) 
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- £ p(«(t+l)/B^(t),l(t))P^(t)  (7.29) 

J«4af  (7,2S)«  (7.2i)»  and  (7.27)  one  obtains  an  nf^ainn 
ali^lar  to  (7.21)/  i.s.* 

Pj(t+1)  - P(H^(t+l)/B(t*fl)) 


M 

■ E MHj/H,^)P(H^/H,5)  MH,^)P,^(t)  (7.30) 

k-1 


^ara  P(B^/Rj^) 


MHj/H^) 

KV 


P<B.(t+l)/H^(t),*(t)) 

transition  probability 

conditional  lilcalihood  ratio  dafined 
in  (7.26) 

likelihood  ratio  defined  in  (7.27) 


The  difference  of  (7.21)  and  (7.30)  is  that  with  limited 
memory/  we  are  interested  in  P(H^  (t-fl)/Z(t+l) ) and  not  in 
P(B^  (t-fl)/Z(t'i-l) ) . This  also  limits  the  number  of  filters 
to  the  number  of  critical  issue  of  this  ap> 

proach  is  the  selection  of  the  transition  probability 
P(H^/Hj^).  In  practical  problems,  it  may  be  selected 
s priori  with  engineering  intuition  and  physical  reasons. 
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(3)  Aging  Filter  Approach 

Nhon  the  system  dynasdcs  are  uncertain  and  changing 
vith  timer  the  aging  filter  {36 » 37]  is  often  used  to 
place  eseponentially  higher  weighting  to  the  more  recent 
raeasureronts . Its  extension  to  the  NMSA  case  (e.g.r  In 
the  probability  computation)  is  not  available.  Prelimi- 
nary results  are  discussed  in  (38] . 

(4)  Others 

There  exist  many  methods  that  can  be  applied  to 
open  up  the  bandwidth  of  each  Kalman  filter  and  to  pre- 
vent the  a posteriori  hypothesis  probability  from  locking 
on  zero  (or  unity) . The  method  used  in  the  previous 
section r i.e.r  increase  process  noise  and  bound  the  prob- 
ability, is  Indeed  just  one  of  them. 

A useful  study  would  be  to  compare  the  above  approaches 
by  applying  them  to  a significant  physical  problem,  such  as  the 
Re-entry  Vehicle  Tracking  problem. 
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7.3  Furthar  Problem  Areaa 


In  thla  aubsaction,  va  conclude  by  suggesting  tha  follow- 
ing further  problem  areas. 

(1)  Prom  section  3.7.  It  was  foxmd  that  8«ae  fundamen- 

• I 

tal  issues  of  MMEA  pertaining  to  its  convergence  and 
identif lability  still  require  rigorous  investigation. 

(2)  It  is  demonstrated  in  section  6 that  a known  input 
may  be  required  in  some  situations  to  help  identify  time- 
varying  parameters.  The  problem  of  optimal  signal  design 
in  using  MMEA  for  system  identification  is  still  an  open 
issue. 

(3)  Further  studies  are  required  to  extend  MMEA  to  time- 
varying  parameters.  The  optimum  MMEA  for  a special  class 
of  time-varying  parameters  and  several  suboptlmal  ap- 
proaches are  discussed  in  section  7.2.  The  extension  of 
MMEA  to  other  types  of  parameter  variation  is  needed.  A 
Comparative  study  of  the  suboptlmal  approaches  is  an 
interesting  further  topic. 
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APPSNDIX  A 

THB  DISCRBTB  KAJLNAN  AND  BXT8NDB0 
KALMAN  PZLTBR  ALGORITHMS 

Zn  this  sppsttdiXf  VS  stats  ths  discrets  Kalnan  flltsr 
algoritha  and  its  first  order  extension  (the  extended  Kaloan 
filter)  to  the  nonlinear  oase. 

A.l  The  Discrete  Kalaan  Filter  Algorithm 

Consider  the  discrete  system  represented  by 

x(t+l)  - A x(t)  + B u(t)  + L £(t)  (A.l) 

%fith  sMasurwoent  equation  represented  by 

8(t+l)  -Cx,(t+1)  + e(t+l)  (A. 2) 

where  x,  and  £ are  state,  control,  and  measurement  vectors, 
respectively,  ^(t)  and  9(t)  are  white  Gaussian  noise  sequences 
with  zero  mean  and  covariances  H(t)  and  6,  respectively.  The 
matrices.  A,  B,  L,  and  C may  be  time-varying  although  not  expli- 
citly shown.  The  discrete  Kalman  filter  algorithm  is  stated 
below. 

Predict  Cycle 

£(t+l/t)  - A £(t/t)  + B u(t)  (A. 3) 

E(t+l/t)  - A E{t/t)A’’  + L 5(t)L'^  (A. 4) 
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Update  Cycle 


x(t+l/t+l)  - x(t+l/t)  + W(t+1) (z(t+l)  -Cx(t+l/t))  (A. 5) 


W(t-H)  - E(t+l/t)C^[C  £(t-H/t)C^  + e(t+l)l“^ 


(A. 6) 


I(t+l/t+l)  « (I  - W<t+l)Cl  Z(t+l/t) 


(A. 7) 


x(t/t)  » E(x(t)/Z(t)) 


(A. 8) 


x(t+l/t)  » E(x(t+1)/Z(t) ) 


(A. 9) 


E (t/t)  » cov(x(t)  ;x(t)/Z(t) ) 


(A. 10) 


E(t+l/t)  » cov(x(t+l),  x(t+l)/Z(t)) 


Z(t)  « the  set  of  all  past  measurements 


(A. 11) 


{u(0),  u(l)  , . . . ,u(t-l) ,z (1) , . . . ,z  (t)  } (A. 12) 


The  initial  estimate  x(0/0)  is  assumed  to  be  Gaussian  with  mean 


x{C)  and  covariance  E(0/0). 


A. 2 The  Discrete  Extended  Kalman  Filter  Algorithm 

i 

Consider  a nonlinear  system  represented  by 


X(t+1)  = f(x(t))  + B U(t);+  L £(t) 


(A. 13) 


with  measurement  equation  represented  by 
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z(t+l)  - h(x(t+D)  + 0(t+i) 


{A. 14) 


where  all  the  matrices  and  vectors  are  the  same  as  previously 
defined  except  that  f ( ) and  h(  ) now  represent  nonlinear  system 
and  measurement  eguatior  'espectively.  The  extended  Kalman 
filter  is  derived  by  expending  £(  ) and  h(  ) in  using  the  Taylor 
series  expansion  up  to  first  order  term.  Let 


P a Jacobian  matrix  of  £ ( ) 


3f  (x(t)) 


x(t)  - x(t/t) 


(A. 15) 


= Jacobian  matrix  of  h(  ) 


3 h(x(t^l)) 

3 xTt-H)  X (t+1)  « X (t+l/t) 


(A. 16) 


The  discrete  extended  Kalman  filter  algorithm  is  stated  below. 


Predict  Cycle 


X (t+l/t)  = f(ic(t/t))  + B u(t) 


E (t+l/t)  » F Z(t/t)F^  + L E(t)L’’ 


(A.17) 


(A. 18) 


Update  Cycle 


fi(t+l/t+l)  * ic(t-^l/t)  + W(t+1)  (s(t+l)  - h(x  (t+l/t)))  (A. 19) 


y(-.n)  « 1 (t+l/t) h'^(h  z(t+i/t)H''’  + e(t+i)i“^  (a.2o) 
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APPENDIX  B 


DISCRETE  LINEAR  SMOOTHING  ALGORITHMS 
.The  system  and  measurement  equations  are  re-stated 

below. 

x(t+l)  - A(t+l,t)X(t)  B U(t)  + L ^(t)  (B.l) 

z(t+l)  - c(t+l)x(t+l)  + e(t+l)  (B.2) 

All  the  definitions  and  statistical  properties  defined  in  the 
Appendix  A still  apply.  Notice  that  the  time-varying  property 
of  A(r+l/t)  and  C(t+1)  is  now  explicitly  shown.  We  still  use  the 
following  definition  for  state  estimate  and  covariance 

,^(T/t)  » Elx(x)/2(t)]  (B.3) 

E(T/t)  « covtx(T)  ;x{T)/Z(t)  ] (B.4) 

Three  kinds  of  smoothing  are  considered.  They  are  de- 
fined below. 

(1)  Fixed-interval  smoothing:  given  Z(T), 
obtain  x(t/T)  and  E(t/T)  for  all  t<T. 

(2)  Fixed-point  smoothing:  given  x, 

obtain  ^ (i/t)  and  £(T/t)  for  all  t>x. 

(3)  Fixed-lag  smoothing:  advance  £(t/t+k) 

and  E(t/t+k)  to  x(t+l/t+l+k)  and 
E(t+l/t+l+k)  where  k is  a positive 
constant. 
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Only  the  algorithms  will  be  stated  here.  Their  deriva- 
tiOM  may  be  found  in  many  references,  e.g.,  refs  [24-28],  These 
algoritlHM  are  stated  individually  in  the  following  subsections. 
•.1  -interval  Smoothing  Algorithm 

In  order  to  use  the  fixed-interval  smoothing  algorithm, 
the  filtering  results  must  be  first  made  available. 

State 


ic(t/T)  » x(t/t)  + G(t)  [«(t+l/T)  - «(t+l/t)) 


(B.5) 


Gain 


G(t)  » Z(t/t)A‘^(t+l,t)5:"^(t-H/t) 


(B.6) 


Covariance 


£(t/T)  = E(t/t)  + G(t) [E(t+1/T)  - E(t+l/t) )G^(t)  (B.7) 


Initial  Conditions 


a (T/T) 


I (T/T) 


(B.8) 


B. 2 Fixed-point  Smoothing  Algorithm 

There  are  several  equivalent  algorithms  in  this  category. 
Only  one  of  them  is  stated  here.  Similary,  the  filterin>7  results 
axe  needed  for  fixed-point  smoothing. 

State 


x(T/t)  -^(T/t-lj,  1 


(B.9) 
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<l>ln 

D(T/t)  • D(T/t-l)E(t/t)A'*’(t+l,t)Z"^(t-H/t)  (B.IO) 

Covariance 

£(T/t)  » E(T/t“l)  - D(x/t)W(t)C(t)l  (t/t-l)D^^(T/t)  (B.ll) 
Initial  Condition 

x(t/t),  E(t/t),  D(tA-)  - I (B.12) 

where  W(t)  » filtering  gain  defined  in  (A. 6). 

I - identity  matrix. 

B.3  Fixed-lag  Smoothing  Algorithm 

In  order  to  perform  fixed-lag  smoothing,  the  filtering, 
fixed-interval  smoothing,  and  fixed-point  smoothing  results  must 
be  available  to  obtain  initial  conditions. 

State 

g(t+l/t+l+k)  - A(t+l,t)£(t/t+k)  + B u(t) 

+ L l£(t/t+k;  - £(t/t)) 

+ D(t+l/t+l+k)W(t+l+k) (z(ttl+k)  - C (t+l+k)x (t+l+k/t+k) ] (B.13) 

Gain 

D(t+l/t+l+k)  - G"^(t)D(t/t+k)G(t+k)  (B.14) 


I -4 


fi 


i/ 


\ 


¥ 


s2 

k 


I 

i 

] 


'i 


1 

■j 

\ 

i 

I 
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Cov>rlanc^ 


E(t+l/t+l+k)  - E(t+l/t)  - D(t4.1/t+l+k) 

• W ( t-i-k-fk ) C ( t+  1+k ) IE  ( t+l+k/ 1 ) ( t + 1/ 1+ 1+k) 

- G“^  (t)-ll(t/t)  - E(t/t+k)JG"^(t)  (B.15) 

where  W(t)  ■ filtering  gain,  defined  in  (A. 6).  G(t)  » flxed**in- 
terval  smoothing  gain,  defined  in  (B.6). 

Initial  Conditions 

x(t  /t  +k),  E(t  /t  +k),  D(t  /t  +k) 

— oo  — oo  — oo 

These  conditions  are  obtained  from  fixed-point  smoothing. 
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